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ABSTRACT 

We study, by means of adaptive mesh refinement hydro- and magnetohydrodynamical simulations that cover a 
wide range of scales (from kpc to sub-parsec), the dimension of the most dissipative structures and the injection 
scale of the turbulent interstellar gas, which we find to be about 75 pc, in agreement with observations. This 
is however smaller than the average size of superbubbles, but consistent with significant density and pressure 
changes in the ISM, which leads to the break-up of bubbles locally and hence to injection of turbulence. The 
scalings of the structure functions are consistent with log-Poisson statistics of supersonic turbulence where energy 
is dissipated mainly through shocks. Our simulations are different from previous ones by other authors as (i) we 
do not assume an isothermal gas, but have temperature variations of several orders of magnitude and (ii) we have 
no artificial forcing of the fluid with some ad hoc Fourier spectrum, but drive turbulence by stellar explosions at 
the Galactic rate, self-regulated by density and temperature thresholds imposed on the ISM gas. 

Subject headings: Hydrodynamics - ISM: structure - Galaxy: evolution - Galaxy: structure - Galaxy: general 
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1. INTRODUCTION 

Long ago it was recognized that the interstellar gas is a highly 
turbulent and compressible medium (von Weizsacker 1951), al- 
though this fact is still being ignored in many theoretical mod- 
els. Indeed, high resolution observations of the ISM show struc- 
tures on all scales down to the smallest resolvable ones, imply- 
ing a dynamical coupling over a wide range of scales, which is 
a main characteristic of a turbulent flow with Reynolds num- 
bers the order of Re = 10 5 - 10 7 (cf. Elmegreen & Scalo 2004). 
It was not until the last decade that the effects of turbulence 
in the ISM in general (e.g., Vazquez-Semadeni et al. 1995, 
2000; Passot et al. 1995; Korpi et al. 1999; Wada & Nor- 
man 2001; Avillez & Mac Low 2002; Avillez & Breitschwerdt 
2004, 2005 - herafter referred to as AB04, AB05; Joung & Mac 
Low 2006) and in molecular clouds, in particular (e.g., Mac 
Low et al. 1998; Stone et al. 1998; Porter et al. 1998, 1999; 
Padoan & Nordlund 1999; Boldyrev et al. 2002; Padoan et al. 
2004), have been studied with high detail by means of numeri- 
cal simulations. The molecular cloud studies focused on decay- 
ing and continuously driven turbulence by means of large-scale 
incompressible forcing in Fourier space, stirring the gas in the 
entire computational domain simultaneously. From these sim- 
ulations it became clear that (i) isothermal supersonic turbulent 
motions inside molecular clouds decay within a sound cross- 
ing time becoming subsonic whether a magnetic field is present 
or not (Mac Low et al. 1998; Stone et al. 1998; Padoan & 
Nordlund 1999) and (ii) the fractal dimension of the most dis- 
sipative structures in molecular clouds evolves from close to 1 
(filaments) to approximately 2 (shocks) with increasing Mach 
number (Padoan et al. 2004). 

Although there has been a lot of observational, theoretical 
and computational research concerning interstellar turbulence, 
little is known about: (i) the scales at which energy is injected, 
and (ii) the dimension of the most dissipative structures, if a 
wide range of scales (from kpc to sub-parsec) is taken into ac- 
count and the turbulence forcing are point explosions acting in 
small regions of the domain as opposed to incompressible forc- 



ing over the entire domain. 

Here we address these issues by means of 3D adaptive mesh 
refinement (AMR) hydro- (HD) and magnetohydrodynamical 
(MHD) simulations of the SN-driven ISM. The outline of the 
present paper is as follows: Section 2 deals with the physical 
model, setup of the simulations and numerical code used; In 
Section 3 a characterisation of the source and injection scale of 
the interstellar turbulence is made, while in Section 4 the Haus- 
dorff dimension of the most dissipative structures is analysed; 
Section 5 closes the paper with a discussion and conclusions. 

2. MODEL AND SIMULATIONS 

We use previously published (AB04, AB05) 3D AMR HD 
and MHD simulations (denoted by HD04 and MHD05) and 
a new 0.625 pc HD run (denoted by HD06) with self-gravity 
of the SN-driven ISM in a section of the Galaxy centered at 
the Solar circle having an area of 1 kpc 2 and extending up to 
10 kpc on either side of the Galactic midplane, capturing both 
the large scale disk-halo-disk circulation as well as the smallest 
structures in the Galactic disk. All the runs use a 10 pc resolu- 
tion coarse grid and generate subgrids up to a local resolution 
of 1.25 pc (3 levels of refinement; HD04 and MHD05 runs) 
and 0.625 pc (4 levels of refinement; HD06 run) in the disk re- 
gion \z\ < 500 pc (corresponding to effective grids of 800 3 and 
1600 3 cells, respectively) and 2.5 pc (2 levels of refinement) on 
the remaining computational domain. The boundary conditions 
are periodic along the vertical faces and outflow on the top and 
bottom faces of the grid. 

The physical model, which does not include heat conduction 
(because we find that turbulent diffusion is a significantly more 
efficient transport process), includes a gravitational field pro- 
vided by the stars in the disk (Kuijken & Gilmore 1989), local 
self-gravity (excluding the contribution from the newly formed 
stars), radiative cooling assuming collisional ionization equilib- 
rium (following Dalgarno & McCray (1972) and Sutherland & 
Dopita (1993) cooling functions for gas with 10 < T < 10 4 and 
10 4 < T < 10 8 5 K, respectively) and solar abundances (Anders 
& Grevesse 1989), uniform heating due to a UV radiation field 
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normalized to the Galactic value and varying with z and photo- 
electric heating of grains and polycyclic aromatic hydrocarbons 
(Wolfire et al. 1995). The disk gas is driven by supernovae 
types la and II (including types lb, Ic and II) with the observed 
Galactic SN rates (Cappellaro et al. 1999). SNe la are set up 
randomly with an exponential distribution having a scale height 
of 325 pc (Freeman 1987). The number and masses of SNe II 
progenitors are determined according to the local mass distri- 
bution in the highest density (n > 10 cm -3 ) cold (T < 100 K) 
regions of converging flows (V • v < 0; here v*is the gas veloc- 
ity) on the level 1 grid (note that the coarse grid is identified as 
level 0) after synchronization of the time steps of the finer and 
coarser grids. In the flagged regions, with a star forming effi- 
ciency of 10%, the available mass that can be turned into stars is 
> 840 Mq. The masses and number of the new stars are deter- 
mined by an initial mass function (IMF) with lower and upper 
mass cutoffs of 8 and 60 Mq, respectively, in order to guarantee 
that all the available mass goes into high mass (> 8M Q ) stars. 
For the present work low mass stars which should be formed 
coevally, are not modelled, but they should have a minor influ- 
ence (except for the stellar gravitational field, which has already 
been accounted for) on the ISM dynamics. 

The newly formed stars, to which the observed mean random 
velocity of 5 km/s is added at time of formation, are followed 
kinematically until the end of their main sequence life time 
determined from Stothers' (1972) formula. The rate at which 
these stars explode is normalized to the Galactic rate, that is, if 
the SN rate in the disk is larger/smaller by 10% of the observed 
value, SN occurrences are reduced/forced artificially. In gen- 
eral 40-50% of the type II SN progenitors explode in the field, 
while the remaining explode correlated in time and space gen- 
erating large superbubble structures, some of which have sizes 
as large as 500 pc (see AB04, AB05). 

The simulations follow a single fluid with an initial verti- 
cal density distribution that includes the contributions due to 
the molecular, neutral, ionized and hot phases observed in the 
ISM (Ferriere 2001). The mean magnetic field component in 
the MHD05 run is initialized assuming equipartition, while the 
random component, SB, is zero. During the first 20 Myr SB 
builds up and the total field becomes 4.5 /iG (AB05). 

We use the 3D Evora- Vienna Astrophysics Fluid-Parallel 
AMR (EVAF-PAMR) Code originally developed by Avillez 
(2000). It is a Fortran 95 code that solves HD and MHD prob- 
lems with AMR (Pember et al. 1996; Balsara 2001) in a paral- 
lel fashion and uses approximate Riemann solvers for the hydro 
and magnetic components (Collela & Woodward 1984; Dai & 
Woodward 1994, 1998). 

3. THE INJECTION SCALE OF INTERSTELLAR TURBULENCE 

The simulations were run for 400 Myr, a time sufficiently 
long to wipe out any signatures of the initial conditions, which 
are still present at 80 Myr of evolution (see AB04). The ini- 
tially unstable vertical gas distribution reaches a global dynam- 
ical equilibrium at ~ 200 Myr after the set-up of the Galactic 
fountain, with hot gas ascending into the halo where it cools 
and finally returns to the disk. During the dynamical equi- 
librium evolution the Galactic midplane is filled with a thin 
cold layer overlayed by a thick frothy disk composed mainly 
of warm neutral and ionized gas having scale heights of 180 
and 900 pc, resembling closely the Lockman (Lockman et al. 
1986) and Reynolds (1987) layers. In the disk the highest den- 
sity (and lowest temperature, if the gas had enough time to cool) 



gas tends to be confined to sheet-like structures (2D), which are 
formed by SN driven shock compressed layers. In the HD runs 
the T < 10 3 K gas has no preferred morphology, while in the 
MHD05 run it tends to be aligned with the magnetic field. The 
highest densities in the present simulations can be up to 10 4 
cm" 3 as a result of the self-gravity occurring in the cold clouds, 
which are mainly formed in regions where several streams of 
convergent flows meet. The orientation of these streams is ran- 
dom. We note in passing, that if clouds are hit by shocks from 
random directions, turbulence in the interior is generated. Tap- 
ping the turbulent ISM energy reservoir, which is huge, rep- 
resents a neat way of sustaining supersonic turbulence inside 
molecular clouds against efficient dissipation. 

The scale at which energy is transferred to the interstellar 
gas can be determined by using the so-called two-point cor- 
relation function Rjj(l,t) = (u,{x + l,t)uj(x,t)), where u,,Uj are 
the components of the fluctuating velocity field it. The diago- 
nal components of Rjj are even functions of / and can be writ- 
ten in terms of the dimensionless scalar functions f(l,t) and 
g(l,t), l=\\l\\, which satisfy f(0,t) = g(0,t) = 1 and /, g < 1, 
as R\\/u 2 = f{l,t) and Rjj/u 2 = g{l,t) if i = j^l and zero if 
i =/ j with u = (| (m)) 1 / 2 . The characteristic injection scale in the 

flow is given by L\\ = J^°° f(l,t)dl, which in the present simu- 
lations is calculated in a region with a linear size of 500 pc at a 
distance of 250 pc from the edges of the computational domain 
to avoid the periodicity effects of the boundary conditions. The 
top panel of Figure Q] shows the history of L\\ in the last 50 
Myr of evolution of the simulated ISM for 1.25 (in black) and 
0.625 pc (in red) resolutions. The time average of L\\ varies 
between 73 and 75 pc for the HD and MHD cases, respectively, 
and seems to be independent of resolution. In all the cases con- 
sidered here, there is a large scatter of L\\ around its mean as 
a result of the expansion of bubbles and superbubbles in an in- 
homogeneous medium characterized by a large number of SNe 
exploding in the field and the merging of bubbles and super- 
bubbles. This scale corresponds to the one at which the power 
spectrum of the solenoidal component of the velocity field has 
a spectral break, which we interpret as a redistribution of en- 
ergy to both larger and smaller scales (Avillez & Breitschwerdt 
2007, in preparation). 

The ratio between the transverse, L22 = C°° g(l,t)dl, and 
longitudinal correlation lengths (L22/Z.11, which in isotropic 
turbulence is 0.5) varies between 0.2 and 1.3 (bottom panel 
of Figure Q~|l, being its time average, over a period of 50 Myr, 
(L22/L1 1 }, ~ 0.5 (HD04 and HD06 runs) and 0.6 (MHD05 run). 
The former value indicates that in a statistical sense the inter- 
stellar unmagnetized turbulence is roughly isotropic, while the 
20% deviation from 0.5 seen in the MHD05 run suggests that 
(L22/L u ), increases with the magnetic field strength. This is 
also supported by similar increases observed in other MHD 
runs we carried out (0.25, 0.5, 1, 1.5 and 2 fiG), which will 
be described in a forthcoming paper. 

4. INTERMITTENCY AND HAUSDORFF DIMENSIONS 

From the Kolmogorov (1941, herafter K41) model for ho- 
mogeneous and isotropic turbulence one can derive that, in the 
inertial range, the moments of order p, (Svf ), are dictated solely 
by (<5v 3 } via the set of scalings C(p)/C(3)» 

(Svf) oc (5v 3 ><»/« 3 \ (1) 
with p > 1 being an integer. In this expression Svi = v(x+l)- 
v(x) with v(x+l) and v(x) being the velocities along the x- axis at 
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two points separated by a distance rj <C / <C L\ \ , with L\ \ being 
the outer scale and r\ = j/ 3 / 4 e -1 / 4 the Kolmogorov microscale; 
( ) stands for the ensemble average over the probability density 
function of <5v/. Experiments have shown that (see discussions 
in Benzi et al. 1993 and Frisch 1995): (i) equation ([TJ, with 
the same scaling exponents, is valid in a wide range of length 
scales for large as well as small Re even if no inertial range is 
established and (ii) in the inertial range of high Re flows the 
scaling exponents become significantly smaller than the linear 
law p/3 (the K41 scaling) with increasing p. This is interpreted 
as the statistics of small scales (the ones dictating the behaviour 
of higher-order moments) becoming increasingly non-gaussian. 

In the unmagnetized and magnetized simulated interstellar 
turbulence driven by point explosions the moments (Svf ) are re- 
lated to (#V; ) (see Figure |2l which shows the moments and the 
slopes C(p)/C(3) °f their best fits for the HD06 run), support- 
ing the validity of equation (Q} in such cases of compressible 
forcing. 

The left panel of Figure [3] compares theoretical, experimen- 
tal and simulated (black and red triangles refer to the HD04 
and HD06 runs, respectively, while squares refer to the MHD05 
run) results on the variation of C(/>)/C(3) with p. The experi- 
mental data is taken from Benzi et al. (1993) and the theoreti- 
cal predictions correspond to the K41 and She-Leveque (1994; 
SL94) models of incompressible turbulence and the Burgers- 
Kolmogorov model (Boldyrev 2002; BK02) for supersonic tur- 
bulence. The K41 model considers a log-normal statistics for 
the transfer of energy from large to small scales and has no cor- 
rections for intermittency. The SL94 model uses a log-Poisson 
statistics to describe intermittency and assumes that energy dis- 
sipation occurs through filamentary structures (Hausdorff di- 
mension D = 1) resulting from vortex stretching. The BK02 
and SL94 models differ in the geometry of the most intermit- 
tent structures, which for supersonic turbulence are assumed to 
be shocks (Hausdorff dimension D = 2). We note that the ge- 
ometry of those scales does not have to represent dissipative 
structures (as in incompressible turbulence) nor do they have to 
be planar. 2D in the Hausdorff sense could represent various 
structures, sheets being the simplest. This is further supported 
by the fact that Burger's turbulence produces discontinuities, 
but not shocks, in the sense that material cannot flow through 
them - they just sweep up material. 

The scalings C(p)/C(3) can be written as (Dubrulle 1994; see 
discussion in Politano & Pouquet 1995) 

^ = £(1 -*) + (3 -Z>)(l (2) 
C(3) g \ ' 

where g is linked to the model of nonlinear energy transfer via 
the basic scaling v; ~ l l l g , x is related to the energy transfer time 
at the smallest scale t\ ~ l x , with x > and (3=1 —x/(3 — D). 
In these expressions D is the fractal dimension of the most in- 
termittent structures. In the present simulations turbulence is 
mostly solenoidal in the inertial range and, therefore, g = 3 and 
x = 2/3. Thus, the Hausdorff dimension D of the most inter- 
mittent structures derived from (f2]l and using the ((p)/((3) val- 
ues shown in the left panel of Figure [3] varies between 1.9 and 
2.2 (D= 1.98-2.02, 1.9-2.03 and 1.91-2.2 in the HD06, 
HD04 and MHD05 runs, respectively; right panel of Figure[3]). 
These values suggest that the energy injected by supernovae 



into interstellar turbulence is dissipated preferentially through 
2D structures that can be identified as shock surfaces in the HD 
cases and also as current sheets in the MHD case (see Miiller 
& Biskamp 2000). The discrepancy observed in the MHD 
run, although within the errors of numerical noise, indicates 
a tendency toward filamentary dissipative structures due to the 
anisotropy induced by the magnetic field. We postpone a more 
detailed discussion on this issue to a forthcoming paper. 

5. DISCUSSION AND CONCLUSIONS 

Contrary to previous work of other authors, the present runs 
incorporate (i) the full extent disk-halo-disk circulation, (ii) 
do not assume isothermality (temperature varies within sev- 
eral orders of magnitude), (iii) include magnetic fields and self- 
gravity, and, (iv) cover a wide range (from kpc to sub-parsec) 
of scales, providing resolution independent information on the 
injection scale, extended self-similarity and dissipative struc- 
tures. The mean integral scale at which energy is injected into 
the interstellar medium is of the order of 75 pc. This scale is sig- 
nificantly lower than the often quoted average size of evolved 
SNRs or superbubbles, because energy is injected already at the 
timescale of break-up of the remnants due to density and pres- 
sure gradients in the inhomogeneous ambient medium, which 
occurs well before they stall. The similarity between the HD 
and MHD results indicates that magnetic pressure and tension 
forces cannot prevent break-out of the hot gas from the bub- 
bles formed by SN explosions, as long as L < /3p£, where L 
and £ are the scale lengths of thermal and magnetic pressures 
(including tension forces) and ftp = 4wP/B 2 is the plasma beta. 
However, the magnetic field introduces a degree of anisotropy 
which is reflected in the 20% deviation from 0.5 seen in the time 
averaged ratio LnjL\\ of the magnetized ISM. Kaplan (1958) 
already suggested an injection scale of 80 pc based on the anal- 
ysis of second order structure functions of interstellar matter. 
The dissipation of energy in our interstellar turbulence simula- 
tions proceeds through shocks and numerical viscosity, as we 
cannot resolve the physical viscosity scales. The variation of 
C(f)/C(3) with P i s most consistent with a log-Poisson model 
for the scales at which intermittency becomes important. 

Recently Joung & Mac Low (2006) found similar scalings 
in a SN-driven ISM HD simulation, although they do not con- 
sider self-gravity, find a hot (T> 10 5 5 K) gas occupation frac- 
tion fh > 40%, which is at odds with observations and other 
models, and their resolution ( 1 .95 pc) in the disk is a factor of 3 
lower than that in the HD06 run (0.625 pc). Their high value of 
fh implies that the amount of Ovi in the ISM is a least a factor 2 
larger than that observed with Copernicus and FUSE. Accord- 
ing to our simulations, reliable comparisons of global proper- 
ties with observations (e.g., occupation fractions) can only be 
made after the system reaches a dynamical equilibrium (f > 150 
Myr), implying a relaxation time much larger than a few million 
years (see discussion in AB04). 

We thank the referee and the editor, John Scalo, for detailed 
reports that allowed us to significantly improve the paper. This 
research is supported by the Portuguese Science Fountation 
(FCT) through project POCTI/FIS/58352/2004. 
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FIG. 1 . — History of the characteristic size (given by L\ \) of the larger eddies (top panel) and of the ratio L,2i/L\ \ (bottom panel) for the MHD (solid line) and HD 
(dashed lines) runs for 1.25 pc (black) and 0.625 pc (red) resolutions. 
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FIG. 2. — Log-log plot of the velocity correlation moments (Svf) as function of (5v? ) and best fits (dashed lines) for p = 2,4, 10. The data refer to the 0.625 pc 
resolution run. 
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FIG. 3. — Left: Exponent £(p)/C(3) for the structure function versus order p. The full line corresponds to £(/>) = p/3 (K41 theory), bullets correspond to data of 
Benzi et al. (1993), red and blue lines refer to the She-Leveque (1994) and Burgers-Kolmogorov (Boldyrev 2002) models, respectively. Triangles represent data 
from HD simulations with 1.25 pc (black) and 0.625 pc (red) resolutions and the green squares refer to the MHD run with 1.25 pc resolution. Fractal dimension D 
of the most dissipative structures (derived from eq.[2) as a function of the exponent f (p)/C(3), with p = 2,4, 10, shown in the left panel for the three simulations. 



